Strategic Analysis of Micromobility Dock Deployment

In [1]:

#install.packages(c("osmdata", "sf", "dplyr", "ggplot2", "rnaturalearth"))
library(osmdata)

Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(sf)

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE

library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(ggplot2)
library(rnaturalearth)
library(osmextract)

Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
Check the package website, https://docs.ropensci.org/osmextract/, for more details.

library(tigris)

To enable caching of data, set `options(tigris_use_cache = TRUE)`
in your R script or .Rprofile.

library(mapview)
library(lehdr)
library(viridis)

Loading required package: viridisLite

library(viridisLite)
library(spatstat)

Loading required package: spatstat.data

Loading required package: spatstat.univar

spatstat.univar 3.1-5

Loading required package: spatstat.geom

spatstat.geom 3.6-1

Loading required package: spatstat.random

spatstat.random 3.4-3

Loading required package: spatstat.explore

Loading required package: nlme

Attaching package: 'nlme'

The following object is masked from 'package:dplyr':

    collapse

spatstat.explore 3.6-0

Loading required package: spatstat.model

Loading required package: rpart

spatstat.model 3.5-0

Loading required package: spatstat.linnet

spatstat.linnet 3.4-0

spatstat 3.5-0 
For an introduction to spatstat, type 'beginner' 

library(units)

udunits database from /usr/share/xml/udunits/udunits2.xml

library(spdep)

Loading required package: spData

To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

library(terra)

terra 1.8.86

Attaching package: 'terra'

The following objects are masked from 'package:spatstat.geom':

    area, delaunay, is.empty, rescale, rotate, shift, where.max,
    where.min

The following object is masked from 'package:tigris':

    blocks

library(raster)

Loading required package: sp

Attaching package: 'raster'

The following object is masked from 'package:nlme':

    getData

The following object is masked from 'package:dplyr':

    select

library(tictoc)

Attaching package: 'tictoc'

The following object is masked from 'package:raster':

    shift

The following objects are masked from 'package:terra':

    shift, size

The following object is masked from 'package:spatstat.geom':

    shift

library(tidycensus)
library(leaflet)

Plotting the Docks

Washington, DC

In [2]:

dc_bb <- c(xmin = -77.45, ymin = 38.70, xmax = -76.80, ymax = 39.10)
docks_raw_dc <- opq(bbox = dc_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Capital Bikeshare") |>
  osmdata_sf()

docks_dc <- docks_raw_dc$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))

In [3]:

hull_dc <- docks_dc |> 
    st_union() |> 
    st_convex_hull()
hull_buffered_dc <- st_buffer(hull_dc, dist = 2500)
hull_buffered_wgs84_dc <- st_transform(hull_buffered_dc, 4326)
docks_dc <- st_intersection(
    st_transform(docks_dc, 4326), 
    hull_buffered_wgs84_dc)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

In [4]:

mapview(hull_buffered_wgs84_dc, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_dc, col.regions = "red")

New York, NY

In [5]:

nyc_bb <- c(xmin = -74.15, ymin = 40.55, xmax = -73.70, ymax = 40.95)
docks_raw_nyc <- opq(bbox = nyc_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Citi Bike") |>
  osmdata_sf()

docks_nyc <- docks_raw_nyc$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))

In [6]:

hull_nyc <- docks_nyc |> 
  st_union() |> 
  st_convex_hull()
hull_buffered_nyc <- st_buffer(hull_nyc, dist = 500)
hull_buffered_wgs84_nyc <- st_transform(hull_buffered_nyc, 4326)
docks_nyc <- st_intersection(
  st_transform(docks_nyc, 4326),
  hull_buffered_wgs84_nyc
)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

In [7]:

mapview(hull_buffered_wgs84_nyc, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_nyc, col.regions = "red")

Boston, MA

In [8]:

bos_bb <- c(-71.20, 42.20, -70.90, 42.45)
docks_raw_bos <- opq(bbox = bos_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Bluebikes") |>
  osmdata_sf()

docks_bos <- docks_raw_bos$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))

In [9]:

hull_bos <- docks_bos |> 
  st_union() |> 
  st_convex_hull()
hull_buffered_bos <- st_buffer(hull_bos, dist = 2500)
hull_buffered_wgs84_bos <- st_transform(hull_buffered_bos, 4326)
docks_bos <- st_intersection(
  st_transform(docks_bos, 4326),
  hull_buffered_wgs84_bos
)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

In [10]:

mapview(hull_buffered_wgs84_bos, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_bos, col.regions = "red")

San Francisco, CA

In [11]:

sf_bb <- c(-122.55, 37.70, -122.35, 37.85)
docks_raw_sf <- opq(bbox = sf_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Bay Wheels") |>
  osmdata_sf()

docks_sf <- docks_raw_sf$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))

In [12]:

hull_sf <- docks_sf |> 
  st_union() |> 
  st_convex_hull()
hull_buffered_sf <- st_buffer(hull_sf, dist = 500)
hull_buffered_wgs84_sf <- st_transform(hull_buffered_sf, 4326)
docks_sf <- st_intersection(
  st_transform(docks_sf, 4326),
  hull_buffered_wgs84_sf
)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

In [13]:

mapview(hull_buffered_wgs84_sf, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_sf, col.regions = "red")

Chicago, IL

In [14]:

chi_bb <- c(-87.80, 41.75, -87.55, 42.00)
docks_raw_chi <- opq(bbox = chi_bb) |>
  add_osm_feature(key = "amenity", value = "bicycle_rental") |>
  add_osm_feature(key = "network", value = "Divvy") |>
  osmdata_sf()

docks_chi <- docks_raw_chi$osm_points |>
  st_as_sf() |>
  st_transform(3857) %>%
  filter(!is.na(network))

In [15]:

hull_chi <- docks_chi |> 
  st_union() |> 
  st_convex_hull()
hull_buffered_chi <- st_buffer(hull_chi, dist = 2500)
hull_buffered_wgs84_chi <- st_transform(hull_buffered_chi, 4326)
docks_chi <- st_intersection(
  st_transform(docks_chi, 4326),
  hull_buffered_wgs84_chi
)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

In [16]:

mapview(hull_buffered_wgs84_chi, col.regions = "lightblue", alpha.regions = 0.3) +
  mapview(docks_chi, col.regions = "red")

Plotting Other Data

Washington, DC

In [17]:

#census_api_key("5e8593d6e289e6ae04c33a019f4951ca2bdd9523", install = TRUE)
dc_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("DC", "MD", "VA"),
  geometry = TRUE
)

Getting data from the 2018-2022 5-year ACS

Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.

Fetching block group data by state and combining the result.

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |======================================================================| 100%

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |======================================================================| 100%

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |======================================================================| 100%

In [18]:

dc_bg <- st_transform(dc_bg, st_crs(hull_buffered_wgs84_dc))

dc_bg <- dc_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

dc_bg_clipped <- st_intersection(dc_bg, hull_buffered_wgs84_dc)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

dc_bg_clipped <- dc_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

dc_bg_density <- dc_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )

In [19]:

max_density <- 20000
dc_bg_density$pop_density_plot <- pmin(dc_bg_density$pop_density, max_density)

In [20]:

dc_bbox <- st_bbox(hull_buffered_wgs84_dc)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_dc <- opq(dc_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_dc))

In [21]:

pal <- leaflet::colorNumeric(viridis(20), domain = dc_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = dc_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_dc,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_dc,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )

New York, NY

In [22]:

nyc_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("NY", "NJ"),
  geometry = TRUE
)

Getting data from the 2018-2022 5-year ACS

Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.

Fetching block group data by state and combining the result.

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |======================================================================| 100%

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |======================================================================| 100%

In [23]:

nyc_bg <- st_transform(nyc_bg, st_crs(hull_buffered_wgs84_nyc))

nyc_bg <- nyc_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

nyc_bg_clipped <- st_intersection(nyc_bg, hull_buffered_wgs84_nyc)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

nyc_bg_clipped <- nyc_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

nyc_bg_density <- nyc_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )

In [24]:

max_density <- 75000
nyc_bg_density$pop_density_plot <- pmin(nyc_bg_density$pop_density, max_density)

In [25]:

nyc_bbox <- st_bbox(hull_buffered_wgs84_nyc)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_nyc <- opq(nyc_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_nyc))

In [26]:

pal <- leaflet::colorNumeric(viridis(20), domain = nyc_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = nyc_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_nyc,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_nyc,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )

Boston, MA

In [27]:

bos_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("MA"),
  geometry = TRUE
)

Getting data from the 2018-2022 5-year ACS

Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |======================================================================| 100%

In [28]:

bos_bg <- st_transform(bos_bg, st_crs(hull_buffered_wgs84_bos))

bos_bg <- bos_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

bos_bg_clipped <- st_intersection(bos_bg, hull_buffered_wgs84_bos)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

bos_bg_clipped <- bos_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

bos_bg_density <- bos_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )

In [29]:

max_density <- 25000
bos_bg_density$pop_density_plot <- pmin(bos_bg_density$pop_density, max_density)

In [30]:

bos_bbox <- st_bbox(hull_buffered_wgs84_bos)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_bos <- opq(bos_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_bos))

In [31]:

pal <- leaflet::colorNumeric(viridis(20), domain = bos_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = bos_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_bos,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_bos,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )

San Francisco, CA

In [32]:

sf_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("CA"),
  geometry = TRUE
)

Getting data from the 2018-2022 5-year ACS

Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |======================================================================| 100%

In [33]:

sf_bg <- st_transform(sf_bg, st_crs(hull_buffered_wgs84_sf))

sf_bg <- sf_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

sf_bg_clipped <- st_intersection(sf_bg, hull_buffered_wgs84_sf)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

sf_bg_clipped <- sf_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

sf_bg_density <- sf_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )

In [34]:

max_density <- 25000
sf_bg_density$pop_density_plot <- pmin(sf_bg_density$pop_density, max_density)

In [35]:

sf_bbox <- st_bbox(hull_buffered_wgs84_sf)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_sf <- opq(sf_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_sf))

In [36]:

pal <- leaflet::colorNumeric(viridis(20), domain = sf_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = sf_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_sf,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_sf,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )

Chicago, IL

In [37]:

chi_bg <- get_acs(
  geography = "block group",
  variables = "B01003_001E",
  year = 2022,
  survey = "acs5",
  state = c("IL"),
  geometry = TRUE
)

Getting data from the 2018-2022 5-year ACS

Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |======================================================================| 100%

In [38]:

chi_bg <- st_transform(chi_bg, st_crs(hull_buffered_wgs84_chi))

chi_bg <- chi_bg %>%
  mutate(area_original = as.numeric(st_area(geometry)))

chi_bg_clipped <- st_intersection(chi_bg, hull_buffered_wgs84_chi)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

chi_bg_clipped <- chi_bg_clipped %>%
  mutate(
    area_clipped = as.numeric(st_area(geometry)),
    pop_adjusted = estimate * (area_clipped / area_original),
    area_km2 = area_clipped / 1e6,
    pop_density = pop_adjusted / area_km2
  )

chi_bg_density <- chi_bg_clipped %>%
  mutate(
    area_km2 = as.numeric(st_area(geometry)) / 1e6,
    pop_density = estimate / area_km2
  )

In [39]:

max_density <- 30000
chi_bg_density$pop_density_plot <- pmin(chi_bg_density$pop_density, max_density)

In [40]:

chi_bbox <- st_bbox(hull_buffered_wgs84_chi)

# bus_stops <- opq(bbox) %>%
#   add_osm_feature(key = "highway", value = "bus_stop") %>%
#   osmdata_sf() %>%
#   .$osm_points %>%
#   st_transform(st_crs(hull_buffered_wgs84_dc))

stations_chi <- opq(chi_bbox) %>%
  add_osm_feature(key = "railway", value = "station") %>%
  osmdata_sf() %>%
  .$osm_points %>%
  st_transform(st_crs(hull_buffered_wgs84_chi))

In [41]:

pal <- leaflet::colorNumeric(viridis(20), domain = chi_bg_density$pop_density_plot)
leaflet() %>%
  leaflet::addProviderTiles("CartoDB.Positron") %>%
  
  # Population density polygons
  leaflet::addPolygons(
    data = chi_bg_density,
    fillColor = ~pal(pop_density_plot),
    fillOpacity = 0.6,
    color = "grey",
    weight = 0.5,
    group = "Population Density",
    popup = ~paste0("Population density: ", round(pop_density_plot))
  ) %>%
  
  # Dock locations
  leaflet::addCircleMarkers(
    data = docks_chi,
    radius = 1,
    color = "blue",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Bikeshare Docks"
  ) %>%
  
  # Train stations
  leaflet::addCircleMarkers(
    data = stations_chi,
    radius = 1,
    color = "red",
    fill = TRUE,
    fillOpacity = 0.5,
    group = "Train Stations"
  ) %>%
  
  # Layer controls
  leaflet::addLayersControl(
    overlayGroups = c("Population Density", "Bikeshare Docks", "Train Stations"),
    options = leaflet::layersControlOptions(collapsed = FALSE)
  )

Univariate Testing

Washington, DC

In [42]:

coords <- st_coordinates(docks_dc)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)

Warning in knn2nb(knn): neighbour object has 7 sub-graphs

lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 35.752, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.8444759043     -0.0012936611      0.0005596277 

In [43]:

dock_values <- inv_dist

localI <- localmoran(dock_values, lw, zero.policy = TRUE)
docks_dc$localI <- localI[,1]

leaflet(docks_dc) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(
    radius = 4,
    color = ~viridis::viridis(20)[cut(localI, breaks = 20)],
    fillOpacity = 0.8,
    popup = ~paste0("Local Moran's I: ", round(localI, 3))
  )

New York, NY

In [44]:

coords <- st_coordinates(docks_nyc)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)

Warning in knn2nb(knn): neighbour object has 2 sub-graphs

lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 48.697, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.6827748494     -0.0004344049      0.0001968368 

Boston, MA

In [45]:

coords <- st_coordinates(docks_bos)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)

lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 23.091, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.842640558      -0.003134796       0.001341649 

San Francisco, CA

In [46]:

coords <- st_coordinates(docks_sf)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)

lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 12.936, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.808442476      -0.009433962       0.003997522 

Chicago, IL

In [47]:

coords <- st_coordinates(docks_chi)

knn <- knearneigh(coords, k = 4)
nb <- knn2nb(knn)

lw <- nb2listw(nb, style = "W")

inv_dist <- sapply(1:nrow(coords), function(i) {
  neighbors <- nb[[i]]
  if(length(neighbors) == 0) return(0)
  1 / mean(sqrt((coords[i,1]-coords[neighbors,1])^2 + (coords[i,2]-coords[neighbors,2])^2))
})

moran_result <- moran.test(inv_dist, lw, zero.policy = TRUE)
moran_result

    Moran I test under randomisation

data:  inv_dist  
weights: lw    

Moran I statistic standard deviate = 43.616, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.9584339511     -0.0011013216      0.0004839866 

Results

Metropolitan Area	Moran’s I Statistic
Washington, DC	\(0.8444759043\)
New York, NY	\(0.6827748494\)
Boston, MA	\(0.842640558\)
San Francisco, CA	\(0.808442476\)
Chicago, IL	\(0.9584339511\)

Hypothesis Test: Population Density

Washington, DC

In [48]:

proj_crs <- "EPSG:32618"
docks_proj <- st_transform(docks_dc, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_dc, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
dc_bg_density_proj <- st_transform(dc_bg_density, 26918)
rtemplate <- rast(
  extent = ext(dc_bg_density_proj),
  resolution = 250,
  crs = crs(dc_bg_density)
)

pop_rast <- rasterize(
  vect(dc_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)

In [49]:

dc_bg_density_proj <- st_transform(dc_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(dc_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(dc_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))

In [50]:

nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)

In [51]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)
summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  774 points
Average intensity 5.06e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 133 vertices
enclosing rectangle: [290970.7, 343621] x [4291156, 4332019] units
                     (52650 x 40860 units)
Window area = 1530020000 square units
Fraction of frame area: 0.711

Dummy quadrature points:
     64 x 64 grid of dummy points, plus 4 corner points
     dummy spacing: 822.6608 x 638.4921 units

Original dummy parameters: =
Planar point pattern:  2974 points
Average intensity 1.94e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 133 vertices
enclosing rectangle: [290970.7, 343621] x [4291156, 4332019] units
                     (52650 x 40860 units)
Window area = 1530020000 square units
Fraction of frame area: 0.711
Quadrature weights:
     (counting weights based on 64 x 64 array of rectangular tiles)
All weights:
    range: [47800, 525000]  total: 1.53e+09
Weights on data points:
    range: [47800, 263000]  total: 1.32e+08
Weights on dummy points:
    range: [47800, 525000]  total: 1.4e+09
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.512591e+01  1.913035e-04 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.512591e+01 4.795913e-02 -1.521990e+01 -1.503191e+01   ***
pop_im       1.913035e-04 5.669327e-06  1.801919e-04  2.024152e-04   ***
                  Zval
(Intercept) -315.39161
pop_im        33.74361

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.512591e+01  1.913035e-04 

Fitted exp(theta):
 (Intercept)       pop_im 
2.697132e-07 1.000191e+00 

anova(m0, m1, test="LR")

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   683.18 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

New York, NY

In [52]:

docks_proj <- st_transform(docks_nyc, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_nyc, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
nyc_bg_density_proj <- st_transform(nyc_bg_density, 26918)
rtemplate <- rast(
  extent = ext(nyc_bg_density_proj),
  resolution = 250,
  crs = crs(nyc_bg_density)
)

pop_rast <- rasterize(
  vect(nyc_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)

In [53]:

nyc_bg_density_proj <- st_transform(nyc_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(nyc_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(nyc_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))

In [54]:

nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)

In [55]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)

Warning: Values of the covariate 'pop_im' were NA or undefined at 4.3% (345 out
of 8083) of the quadrature points. Occurred while executing: ppm.ppp(Q =
dock_ppp, trend = ~pop_im, data = NULL, interaction = NULL)

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  2303 points
Average intensity 5.99e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 137 vertices
enclosing rectangle: [576179.5, 597719.6] x [4495675, 4527093] units
                     (21540 x 31420 units)
Window area = 384647000 square units
Fraction of frame area: 0.568

Dummy quadrature points:
     100 x 100 grid of dummy points, plus 4 corner points
     dummy spacing: 215.4011 x 314.1755 units

Original dummy parameters: =
Planar point pattern:  5780 points
Average intensity 1.5e-05 points per square unit
Window: polygonal boundary
single connected closed polygon with 137 vertices
enclosing rectangle: [576179.5, 597719.6] x [4495675, 4527093] units
                     (21540 x 31420 units)
Window area = 384647000 square units
Fraction of frame area: 0.568
Quadrature weights:
     (counting weights based on 100 x 100 array of rectangular tiles)
All weights:
    range: [8530, 67700]    total: 3.84e+08
Weights on data points:
    range: [11300, 33800]   total: 69200000
Weights on dummy points:
    range: [8530, 67700]    total: 3.15e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.237567e+01  2.121081e-05 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.237567e+01 3.074223e-02 -1.243592e+01 -1.231541e+01   ***
pop_im       2.121081e-05 9.451471e-07  1.935836e-05  2.306327e-05   ***
                  Zval
(Intercept) -402.56236
pop_im        22.44181

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.237567e+01  2.121081e-05 

Fitted exp(theta):
 (Intercept)       pop_im 
4.220043e-06 1.000021e+00 
Problem:
 Values of the covariate 'pop_im' were NA or undefined at 4.3% (345 out of 8083) of the quadrature points 

anova(m0, m1, test="LR")

Warning: Values of the covariate 'pop_im' were NA or undefined at 4.3% (345 out
of 8083) of the quadrature points. Occurred while executing: ppm.ppp(Q =
dock_ppp, trend = ~pop_im, data = NULL, interaction = NULL,

Warning: Models were re-fitted after discarding quadrature points that were
illegal under some of the models

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance  Pr(>Chi)    
1  346                          
2  347  1   436.56 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Boston, MA

In [56]:

docks_proj <- st_transform(docks_bos, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_bos, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
bos_bg_density_proj <- st_transform(bos_bg_density, 26918)
rtemplate <- rast(
  extent = ext(bos_bg_density_proj),
  resolution = 250,
  crs = crs(bos_bg_density)
)

pop_rast <- rasterize(
  vect(bos_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)

In [57]:

bos_bg_density_proj <- st_transform(bos_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(bos_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(bos_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))

In [58]:

nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)

In [59]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)

Warning: Values of the covariate 'pop_im' were NA or undefined at 4.4% (145 out
of 3260) of the quadrature points. Occurred while executing: ppm.ppp(Q =
dock_ppp, trend = ~pop_im, data = NULL, interaction = NULL)

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  320 points
Average intensity 1.17e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 127 vertices
enclosing rectangle: [812090.4, 831648.2] x [4686171, 4706029] units
                     (19560 x 19860 units)
Window area = 272770000 square units
Fraction of frame area: 0.702

Dummy quadrature points:
     64 x 64 grid of dummy points, plus 4 corner points
     dummy spacing: 305.5909 x 310.2755 units

Original dummy parameters: =
Planar point pattern:  2940 points
Average intensity 1.08e-05 points per square unit
Window: polygonal boundary
single connected closed polygon with 127 vertices
enclosing rectangle: [812090.4, 831648.2] x [4686171, 4706029] units
                     (19560 x 19860 units)
Window area = 272770000 square units
Fraction of frame area: 0.702
Quadrature weights:
     (counting weights based on 64 x 64 array of rectangular tiles)
All weights:
    range: [12100, 94800]   total: 2.72e+08
Weights on data points:
    range: [19000, 47400]   total: 13400000
Weights on dummy points:
    range: [12100, 94800]   total: 2.59e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.418988e+01  9.188116e-05 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.418988e+01 8.618496e-02 -1.435880e+01 -1.402096e+01   ***
pop_im       9.188116e-05 8.049611e-06  7.610421e-05  1.076581e-04   ***
                  Zval
(Intercept) -164.64449
pop_im        11.41436

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.418988e+01  9.188116e-05 

Fitted exp(theta):
 (Intercept)       pop_im 
6.877236e-07 1.000092e+00 
Problem:
 Values of the covariate 'pop_im' were NA or undefined at 4.4% (145 out of 3260) of the quadrature points 

anova(m0, m1, test="LR")

Warning: Values of the covariate 'pop_im' were NA or undefined at 4.4% (145 out
of 3260) of the quadrature points. Occurred while executing: ppm.ppp(Q =
dock_ppp, trend = ~pop_im, data = NULL, interaction = NULL,

Warning: Models were re-fitted after discarding quadrature points that were
illegal under some of the models

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance  Pr(>Chi)    
1  146                          
2  147  1    105.4 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

San Francisco, CA

In [60]:

docks_proj <- st_transform(docks_sf, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_sf, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
sf_bg_density_proj <- st_transform(sf_bg_density, 26918)
rtemplate <- rast(
  extent = ext(sf_bg_density_proj),
  resolution = 250,
  crs = crs(sf_bg_density)
)

pop_rast <- rasterize(
  vect(sf_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)

In [61]:

sf_bg_density_proj <- st_transform(sf_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(sf_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(sf_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))

In [62]:

nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)

In [63]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)

Warning: Values of the covariate 'pop_im' were NA or undefined at 0.66% (6 out
of 912) of the quadrature points. Occurred while executing: ppm.ppp(Q =
dock_ppp, trend = ~pop_im, data = NULL, interaction = NULL)

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  107 points
Average intensity 8.44e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 136 vertices
enclosing rectangle: [-3757195, -3744439] x [5412459, 5425830] units
                     (12760 x 13370 units)
Window area = 126726000 square units
Fraction of frame area: 0.743

Dummy quadrature points:
     32 x 32 grid of dummy points, plus 4 corner points
     dummy spacing: 398.6460 x 417.8386 units

Original dummy parameters: =
Planar point pattern:  805 points
Average intensity 6.35e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 136 vertices
enclosing rectangle: [-3757195, -3744439] x [5412459, 5425830] units
                     (12760 x 13370 units)
Window area = 126726000 square units
Fraction of frame area: 0.743
Quadrature weights:
     (counting weights based on 32 x 32 array of rectangular tiles)
All weights:
    range: [7200, 167000]   total: 1.27e+08
Weights on data points:
    range: [41600, 83300]   total: 8120000
Weights on dummy points:
    range: [7200, 167000]   total: 1.19e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.432761e+01  3.904106e-05 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.432761e+01 1.743289e-01 -1.466929e+01 -1.398594e+01   ***
pop_im       3.904106e-05 1.500172e-05  9.638223e-06  6.844389e-05    **
                  Zval
(Intercept) -82.187281
pop_im        2.602438

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.432761e+01  3.904106e-05 

Fitted exp(theta):
 (Intercept)       pop_im 
5.992334e-07 1.000039e+00 
Problem:
 Values of the covariate 'pop_im' were NA or undefined at 0.66% (6 out of 912) of the quadrature points 

anova(m0, m1, test="LR")

Warning: Values of the covariate 'pop_im' were NA or undefined at 0.66% (6 out
of 912) of the quadrature points. Occurred while executing: ppm.ppp(Q =
dock_ppp, trend = ~pop_im, data = NULL, interaction = NULL,

Warning: Models were re-fitted after discarding quadrature points that were
illegal under some of the models

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance Pr(>Chi)  
1    7                       
2    8  1   6.4363  0.01118 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Chicago, IL

In [64]:

docks_proj <- st_transform(docks_chi, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_chi, 32618)
boundary_owin <- as.owin(boundary_proj)

dock_coords <- st_coordinates(docks_proj)

dock_ppp <- ppp(
  x = dock_coords[,1],
  y = dock_coords[,2],
  window = boundary_owin
)
chi_bg_density_proj <- st_transform(chi_bg_density, 26918)
rtemplate <- rast(
  extent = ext(chi_bg_density_proj),
  resolution = 250,
  crs = crs(chi_bg_density)
)

pop_rast <- rasterize(
  vect(chi_bg_density_proj),
  rtemplate,
  field = "pop_density_plot"
)

In [65]:

chi_bg_density_proj <- st_transform(chi_bg_density, crs = proj_crs)

rtemplate <- rast(
  extent = ext(chi_bg_density_proj),
  resolution = 250,
  crs = "EPSG:32618"
)

pop_rast <- rasterize(
  vect(chi_bg_density_proj),
  rtemplate,
  field = "pop_density_plot",
  fun = "mean"
)
pop_rast <- mask(pop_rast, vect(boundary_proj))

In [66]:

nx <- ncol(pop_rast)
ny <- nrow(pop_rast)
r_ext <- ext(pop_rast)
res_x <- res(pop_rast)[1]
res_y <- res(pop_rast)[2]

xcol <- seq(r_ext$xmin + res_x/2, r_ext$xmax - res_x/2, length.out = nx)
yrow <- seq(r_ext$ymin + res_y/2, r_ext$ymax - res_y/2, length.out = ny)

vals <- values(pop_rast, mat = FALSE)
mat <- matrix(vals, nrow = ny, ncol = nx, byrow = TRUE)

mat_oriented <- mat[ny:1, , drop = FALSE]
pop_im <- im(mat_oriented, xcol, yrow)

In [67]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ pop_im)

Warning: Values of the covariate 'pop_im' were NA or undefined at 10% (444 out
of 4391) of the quadrature points. Occurred while executing: ppm.ppp(Q =
dock_ppp, trend = ~pop_im, data = NULL, interaction = NULL)

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ pop_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  909 points
Average intensity 9.85e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 131 vertices
enclosing rectangle: [-566382.1, -541021] x [4685608, 4737799] units
                     (25360 x 52190 units)
Window area = 922506000 square units
Fraction of frame area: 0.697

Dummy quadrature points:
     70 x 70 grid of dummy points, plus 4 corner points
     dummy spacing: 362.3008 x 745.5848 units

Original dummy parameters: =
Planar point pattern:  3482 points
Average intensity 3.77e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 131 vertices
enclosing rectangle: [-566382.1, -541021] x [4685608, 4737799] units
                     (25360 x 52190 units)
Window area = 922506000 square units
Fraction of frame area: 0.697
Quadrature weights:
     (counting weights based on 70 x 70 array of rectangular tiles)
All weights:
    range: [24600, 270000]  total: 9.21e+08
Weights on data points:
    range: [24600, 135000]  total: 1.08e+08
Weights on dummy points:
    range: [24600, 270000]  total: 8.14e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~pop_im
Model depends on external covariate 'pop_im'
Covariates provided:
    pop_im: im

Fitted trend coefficients:
  (Intercept)        pop_im 
-1.419576e+01  9.193927e-05 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.419576e+01 4.680314e-02 -1.428750e+01 -1.410403e+01   ***
pop_im       9.193927e-05 4.699411e-06  8.272859e-05  1.011499e-04   ***
                 Zval
(Intercept) -303.3079
pop_im        19.5640

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)        pop_im 
-1.419576e+01  9.193927e-05 

Fitted exp(theta):
 (Intercept)       pop_im 
6.836891e-07 1.000092e+00 
Problem:
 Values of the covariate 'pop_im' were NA or undefined at 10% (444 out of 4391) of the quadrature points 

anova(m0, m1, test="LR")

Warning: Values of the covariate 'pop_im' were NA or undefined at 10% (444 out
of 4391) of the quadrature points. Occurred while executing: ppm.ppp(Q =
dock_ppp, trend = ~pop_im, data = NULL, interaction = NULL,

Warning: Models were re-fitted after discarding quadrature points that were
illegal under some of the models

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~pop_im     Poisson
  Npar Df Deviance  Pr(>Chi)    
1  445                          
2  446  1   269.89 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Results

City	Docks	Dock Intensity (points/unit squared)	pop_im coefficient	Z-value	Deviance diff
Washington, DC	\(774\)	\(5.06 \times 10^{-7}\)	\(1.913 \times 10^{-4}\)	\(33.74\)	\(683.18\)
New York, NY	\(2,303\)	\(1.5 \times 10^{-5}\)	\(2.121 \times 10^-5\)	\(22.44\)	\(436.56\)
Boston, MA	\(320\)	\(1.08 \times 10^{-5}\)	\(9.189 \times 10^{-5}\)	\(11.41\)	\(105.4\)
San Francisco, CA	\(107\)	\(6.35 \times 10^{-6}\)	\(3.904 \times 10^{-5}\)	\(2.60\)	\(6.436\)
Chicago, IL	\(909\)	\(3.77 \times 10^{-6}\)	\(9.194 \times 10^{-5}\)	\(19.56\)	\(269.89\)

Hypothesis Test: Other Transit options

Washington, DC

In [68]:

stations_proj <- st_transform(stations_dc, 32618)
docks_proj <- st_transform(docks_dc, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_dc, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)

Warning: 7 points were rejected as lying outside the specified window

dist_im <- distmap(station_ppp)

In [69]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   1112.3 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  774 points
Average intensity 5.06e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 133 vertices
enclosing rectangle: [290970.7, 343621] x [4291156, 4332019] units
                     (52650 x 40860 units)
Window area = 1530020000 square units
Fraction of frame area: 0.711

Dummy quadrature points:
     64 x 64 grid of dummy points, plus 4 corner points
     dummy spacing: 822.6608 x 638.4921 units

Original dummy parameters: =
Planar point pattern:  2974 points
Average intensity 1.94e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 133 vertices
enclosing rectangle: [290970.7, 343621] x [4291156, 4332019] units
                     (52650 x 40860 units)
Window area = 1530020000 square units
Fraction of frame area: 0.711
Quadrature weights:
     (counting weights based on 64 x 64 array of rectangular tiles)
All weights:
    range: [47800, 525000]  total: 1.53e+09
Weights on data points:
    range: [47800, 263000]  total: 1.32e+08
Weights on dummy points:
    range: [47800, 525000]  total: 1.4e+09
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-12.611735447  -0.001085935 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -12.611735447 0.0604945559 -12.730302598 -12.493168296   ***
dist_im      -0.001085935 0.0000453318  -0.001174783  -0.000997086   ***
                  Zval
(Intercept) -208.47720
dist_im      -23.95525

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-12.611735447  -0.001085935 

Fitted exp(theta):
 (Intercept)      dist_im 
3.332674e-06 9.989147e-01 

New York, NY

In [70]:

stations_proj <- st_transform(stations_nyc, 32618)
docks_proj <- st_transform(docks_nyc, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_nyc, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)

Warning: 75 points were rejected as lying outside the specified window

dist_im <- distmap(station_ppp)

In [71]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   1174.6 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  2303 points
Average intensity 5.99e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 137 vertices
enclosing rectangle: [576179.5, 597719.6] x [4495675, 4527093] units
                     (21540 x 31420 units)
Window area = 384647000 square units
Fraction of frame area: 0.568

Dummy quadrature points:
     100 x 100 grid of dummy points, plus 4 corner points
     dummy spacing: 215.4011 x 314.1755 units

Original dummy parameters: =
Planar point pattern:  5780 points
Average intensity 1.5e-05 points per square unit
Window: polygonal boundary
single connected closed polygon with 137 vertices
enclosing rectangle: [576179.5, 597719.6] x [4495675, 4527093] units
                     (21540 x 31420 units)
Window area = 384647000 square units
Fraction of frame area: 0.568
Quadrature weights:
     (counting weights based on 100 x 100 array of rectangular tiles)
All weights:
    range: [8530, 67700]    total: 3.84e+08
Weights on data points:
    range: [11300, 33800]   total: 69200000
Weights on dummy points:
    range: [8530, 67700]    total: 3.15e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-11.074435702  -0.001612625 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -11.074435702 3.398891e-02 -11.141052740 -11.007818663   ***
dist_im      -0.001612625 6.225954e-05  -0.001734652  -0.001490599   ***
                  Zval
(Intercept) -325.82498
dist_im      -25.90166

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-11.074435702  -0.001612625 

Fitted exp(theta):
 (Intercept)      dist_im 
1.550364e-05 9.983887e-01 

Boston, MA

In [72]:

stations_proj <- st_transform(stations_bos, 32618)
docks_proj <- st_transform(docks_bos, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_bos, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)

Warning: 25 points were rejected as lying outside the specified window

dist_im <- distmap(station_ppp)

In [73]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   263.95 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  320 points
Average intensity 1.17e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 127 vertices
enclosing rectangle: [812090.4, 831648.2] x [4686171, 4706029] units
                     (19560 x 19860 units)
Window area = 272770000 square units
Fraction of frame area: 0.702

Dummy quadrature points:
     64 x 64 grid of dummy points, plus 4 corner points
     dummy spacing: 305.5909 x 310.2755 units

Original dummy parameters: =
Planar point pattern:  2940 points
Average intensity 1.08e-05 points per square unit
Window: polygonal boundary
single connected closed polygon with 127 vertices
enclosing rectangle: [812090.4, 831648.2] x [4686171, 4706029] units
                     (19560 x 19860 units)
Window area = 272770000 square units
Fraction of frame area: 0.702
Quadrature weights:
     (counting weights based on 64 x 64 array of rectangular tiles)
All weights:
    range: [12100, 94800]   total: 2.72e+08
Weights on data points:
    range: [19000, 47400]   total: 13400000
Weights on dummy points:
    range: [12100, 94800]   total: 2.59e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-12.306690343  -0.001460137 

                 Estimate         S.E.     CI95.lo       CI95.hi Ztest
(Intercept) -12.306690343 0.0944239275 -12.4917578 -12.121622846   ***
dist_im      -0.001460137 0.0001187592  -0.0016929  -0.001227373   ***
                  Zval
(Intercept) -130.33445
dist_im      -12.29494

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-12.306690343  -0.001460137 

Fitted exp(theta):
 (Intercept)      dist_im 
4.521393e-06 9.985409e-01 

San Francisco, CA

In [74]:

stations_proj <- st_transform(stations_sf, 32618)
docks_proj <- st_transform(docks_sf, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_sf, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)

Warning: 14 points were rejected as lying outside the specified window

dist_im <- distmap(station_ppp)

In [75]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance Pr(>Chi)   
1    1                        
2    2  1    9.655 0.001888 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  107 points
Average intensity 8.44e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 136 vertices
enclosing rectangle: [-3757195, -3744439] x [5412459, 5425830] units
                     (12760 x 13370 units)
Window area = 126726000 square units
Fraction of frame area: 0.743

Dummy quadrature points:
     32 x 32 grid of dummy points, plus 4 corner points
     dummy spacing: 398.6460 x 417.8386 units

Original dummy parameters: =
Planar point pattern:  805 points
Average intensity 6.35e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 136 vertices
enclosing rectangle: [-3757195, -3744439] x [5412459, 5425830] units
                     (12760 x 13370 units)
Window area = 126726000 square units
Fraction of frame area: 0.743
Quadrature weights:
     (counting weights based on 32 x 32 array of rectangular tiles)
All weights:
    range: [7200, 167000]   total: 1.27e+08
Weights on data points:
    range: [41600, 83300]   total: 8120000
Weights on dummy points:
    range: [7200, 167000]   total: 1.19e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-1.362773e+01 -2.971307e-04 

                 Estimate         S.E.       CI95.lo       CI95.hi Ztest
(Intercept) -1.362773e+01 0.1466369755 -13.915136434 -1.334033e+01   ***
dist_im     -2.971307e-04 0.0001051373  -0.000503196 -9.106536e-05    **
                  Zval
(Intercept) -92.935177
dist_im      -2.826121

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-1.362773e+01 -2.971307e-04 

Fitted exp(theta):
 (Intercept)      dist_im 
1.206565e-06 9.997029e-01 

Chicago, IL

In [76]:

stations_proj <- st_transform(stations_chi, 32618)
docks_proj <- st_transform(docks_chi, 32618)
boundary_proj <- st_transform(hull_buffered_wgs84_chi, 32618)

boundary_owin <- as.owin(boundary_proj)
dock_coords <- st_coordinates(docks_proj)
dock_ppp <- ppp(x = dock_coords[,1],
                y = dock_coords[,2],
                window = boundary_owin)

station_coords <- st_coordinates(stations_proj)
station_ppp <- ppp(x = station_coords[,1],
                   y = station_coords[,2],
                   window = boundary_owin)

Warning: 20 points were rejected as lying outside the specified window

dist_im <- distmap(station_ppp)

In [77]:

m0 <- ppm(dock_ppp ~ 1)
m1 <- ppm(dock_ppp ~ dist_im)

anova(m0, m1, test = "LR")

Analysis of Deviance Table

Model 1: ~1      Poisson
Model 2: ~dist_im    Poisson
  Npar Df Deviance  Pr(>Chi)    
1    1                          
2    2  1   386.02 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m1)

Point process model
Fitted to data: dock_ppp
Fitting method: maximum likelihood (Berman-Turner approximation)
Model was fitted using glm()
Algorithm converged
Call:
ppm.formula(Q = dock_ppp ~ dist_im)
Edge correction: "border"
    [border correction distance r = 0 ]
--------------------------------------------------------------------------------
Quadrature scheme (Berman-Turner) = data + dummy + weights

Data pattern:
Planar point pattern:  909 points
Average intensity 9.85e-07 points per square unit
Window: polygonal boundary
single connected closed polygon with 131 vertices
enclosing rectangle: [-566382.1, -541021] x [4685608, 4737799] units
                     (25360 x 52190 units)
Window area = 922506000 square units
Fraction of frame area: 0.697

Dummy quadrature points:
     70 x 70 grid of dummy points, plus 4 corner points
     dummy spacing: 362.3008 x 745.5848 units

Original dummy parameters: =
Planar point pattern:  3482 points
Average intensity 3.77e-06 points per square unit
Window: polygonal boundary
single connected closed polygon with 131 vertices
enclosing rectangle: [-566382.1, -541021] x [4685608, 4737799] units
                     (25360 x 52190 units)
Window area = 922506000 square units
Fraction of frame area: 0.697
Quadrature weights:
     (counting weights based on 70 x 70 array of rectangular tiles)
All weights:
    range: [24600, 270000]  total: 9.21e+08
Weights on data points:
    range: [24600, 135000]  total: 1.08e+08
Weights on dummy points:
    range: [24600, 270000]  total: 8.14e+08
--------------------------------------------------------------------------------
FITTED :

Nonstationary Poisson process

---- Intensity: ----

Log intensity: ~dist_im
Model depends on external covariate 'dist_im'
Covariates provided:
    dist_im: im

Fitted trend coefficients:
  (Intercept)       dist_im 
-1.285494e+01 -8.244326e-04 

                 Estimate         S.E.       CI95.lo      CI95.hi Ztest
(Intercept) -1.285494e+01 5.571363e-02 -1.296414e+01 -12.74574347   ***
dist_im     -8.244326e-04 4.844097e-05 -9.193751e-04  -0.00072949   ***
                  Zval
(Intercept) -230.73239
dist_im      -17.01933

----------- gory details -----

Fitted regular parameters (theta):
  (Intercept)       dist_im 
-1.285494e+01 -8.244326e-04 

Fitted exp(theta):
 (Intercept)      dist_im 
2.613187e-06 9.991759e-01 

Results

City	Docks	Dock Intensity (points/unit squared)	dist_im coefficient	Z-value	Deviance diff
Washington, DC	\(774\)	\(5.06 \times 10^{-7}\)	\(-1.09 \times 10^{-3}\)	\(-23.96\)	\(1112.3\)
New York, NY	\(2,303\)	\(1.5 \times 10^{-5}\)	\(-1.61 \times 10^{-3}\)	\(-25.90\)	\(1174.6\)
Boston, MA	\(320\)	\(1.08 \times 10^{-5}\)	\(-1.46 \times 10^{-3}\)	\(-12.29\)	\(263.95\)
San Francisco, CA	\(107\)	\(8.44 \times 10^{-7}\)	\(-2.97 \times 10^{-04}\)	\(-2.83\)	\(9.655\)
Chicago, IL	\(909\)	\(3.77 \times 10^{-6}\)	\(-8.24 \times 10^{-4}\)	\(-17.02\)	\(386.02\)